//
// SofferBoundProkudin.cxx
//
// Created by TB on 11/23/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <algorithm>
#include <cmath>
#include <cstdlib> // For abs(int)
#include <stdexcept>

#include <boost/foreach.hpp>

#include "tppmc/Particle.h"
#include "tppmc/SofferBoundProkudin.h"
#include "tppmc/PartonContent.h"

// Access to external FORTRAN PARPOLT subroutine.
extern "C" {
   void parpolt_(int *iset, // Order (LO or NLO)
                 double *x, // x
                 double *q2, // scale
                 double *u, // u quark
                 double *d, // d
                 double *ub, // anti-u
                 double *db, // anti-d
                 double *st); // strange/anti-strange
}

namespace tppmc {

   //
   // class SofferBoundProkudin
   //
   
   SofferBoundProkudin::SofferBoundProkudin(int order)
   : mOrder(order) {
   }
   
   SofferBoundProkudin::~SofferBoundProkudin() {
   }

   void SofferBoundProkudin::Evaluate(const pythpart& parton, double Q2,
                                      PartonContent& partcontent) const {
      double  u = 0., d = 0., ub = 0., db = 0., st = 0.;
      double x = parton.X();
      // partol_ needs a non-const integer argument for order
      int order = mOrder;
      // The subroutine returns x * transversity for each flavour
      parpolt_(&order, &x, &Q2, &u, &d, &ub, &db, &st);
      // Initialise to zeros:
      partcontent = PartonContent(0.);
      // Assign parton densities, dividing out the factor of x:
      partcontent[LHAPDF::UP]          =  u / x; // 
      partcontent[LHAPDF::DOWN]         =  (-1.) * d / x; // 
      partcontent[LHAPDF::UBAR]      =  ub / x; //
      partcontent[LHAPDF::DBAR]    =  ub / x; //
      partcontent[LHAPDF::STRANGE]      =  st / x; //
      partcontent[LHAPDF::SBAR] =  st / x; //
   }

   bool SofferBoundProkudin::Supports(const pythpart& parton) const {
      // Light quarks only
      return ::abs(parton.Pdg()) < 4;
   }
}
